### The packages (repetition) ###
   library(here)
   library(sp)
   library(sf)
   library(tmap)
   library(dismo)
   library(raster)
   library(GISTools)
   library(rgdal)
   library(maptools)
   library(viridis)
   library(rgeos)
   library(rgl)
 # library(rasterVis) 
   library(unmarked)
   library(spatstat)

## check that the package terra is detached!
# detach(package:terra)

### The workspace  (repetition) ###
getwd() # you can also use the package 'here()'

# my data are outside my course-folder, therefore I have to do it the old-fashioned way.
root_wd <- "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub"
# relative to work-wd
maps_wd <- paste(root_wd,"/","_DataBorneo/BaseMaps",sep='')
recs_wd <- paste(root_wd,"/","_DataBorneo/BaseRecords",sep='')

# the output folder should have been created by you during Tutorial 2 'R goes spatial'.
# It should contain the hillshade.asc
output_wd <- paste(root_wd,"/","output",sep='')

setwd(output_wd) ## set to the OUTPUT folder!
getwd()

### Plot to get an impression ###

plot(ras42, col=grey.colors(20))
plot(PA_shp,border='green', lwd=1.8, add=T)
plot(sp_recs_sf$geometry, pch= '*',cex=1,col='deeppink',add=T)
text(112, 5, 'Starting with SDMs', cex=2, col= 'red')
plot(River_shp[,3], col='dodgerblue4', add=T)

# workaround for slow computers. First, aggregate the 1 km² resolution into 50*50 km  
agg_pred <- aggregate(x=predictors,fact=50,FUN=mean)
plot(agg_pred)

head(Borneo_shp_sf)

plot(Borneo_shp_sf[,3]) # column 3 is NAME_0

### *Rasterize* the shapefiles
## command rasterize() 
## field = 1 means all raster cells get value = 1, other cell values are set to 0
## background sets all other cells to a chosen number


Mal_ras <- rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_0 %in%
                       c("Brunei","Malaysia"),], y=ras24, field=1,
                       background=0.1)

plot(Mal_ras)

# almost same as:
# NotInd_ras <- rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_0 !=
#                           "Indonesia",], y=ras24, field=1, background=0.01)

Sab_ras <- rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_1 == 'Sabah',],
                     y=ras24, field=1, background=0.1)

PA_ras <- rasterize(x=PA_shp, y=ras24, field=1, background=0.1)


### Crosscheck with plot
par(mfrow=c(1,2))
plot(Sab_ras)
plot(PA_ras) #n.b.: only 1 and small values, no NA

par(mfrow=c(1,1))

ras24_300 <- ras24 <= 300
plot(ras24_300)

bias_tmp <- ras24_300 + Sab_ras + Mal_ras + PA_ras
bias_tmp

plot(bias_tmp)

# get the maximum value in the layer, standardize it and round to two decimals
maxval <- max(values(bias_tmp),na.rm=T)
bias_tmp2 <- bias_tmp/maxval
bias_tmp3 <- round(bias_tmp2, digits=2)

## same as in one go:   
bias1 <- round(bias_tmp/max(values(bias_tmp),na.rm=T),digits=2)

table(values(bias1)) 

# table(bias1@data@values) # is doing the same

### in case of having 0 somewhere:
bias1[values(bias1 == 0)] <- 0.01 # because MaxEnt
                                  # does not take 0!
table(values(bias1))

compareRaster(ras24,bias1) # the same?

# Are they really the same???
length(which(!is.na(values(bias1))))

length(which(!is.na(values(ras24))))

### Plot the bias file
plot(bias1, col = viridis(7)) #or col = grey.colors(7)

# Have a look at the MaxEnt layers produced
# MaxEnt saves the mean of x repetitions in _avg.asc

# The script now only works with the folder structure
# MaxEntRes/nobias
# MaxEntres/withbias

# new argument 'recursive' means, that also subfolders are checked!
infiles <- list.files(path=paste(output_wd,'/MaxEntRes',
                                 sep=''),pattern='_avg.asc$',
                      full.names=TRUE,recursive=TRUE )
infiles

# example
#[1] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/river_dummy_avg.asc"
#[2] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/river_dummy_FutureMaps_avg.asc"
#[3] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/river_dummy_avg.asc"
#[4] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/river_dummy_FutureMaps_avg.asc"

## IF you have a slow computer, 
## please use the workaround with function aggregate() from above:
# me_stack <- aggregate(stack(infiles[c(1:4)]),fact= 50, FUN = mean )

me_stack <- stack(infiles[c(1:4)])

# name sequence according to infiles list above
names(me_stack) <- c('curr_noBias','fut_noBias',
                     'curr_withBias','fut_withBias')

#par(mar=c(4,4,3,3), oma=c(2,2,2,2))

plot(me_stack,col=viridis(100)) # note: I might use a different example than you did in MaxEnt

boxplot(me_stack,layers = c(1,3,2,4),notch=T,outline=F)

# Use maxentResults.csv to access threshold value
# column: 10th percentile training presence logistic threshold

me_res <- list.files(path=paste(output_wd,'/MaxEntRes',sep=''),
                     pattern='maxentResults.csv',
                     full.names=TRUE,recursive=TRUE )
me_res

# example
#[1] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/maxentResults.csv"
#[2] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/maxentResults.csv"

store_res <- lapply(me_res,FUN=read.csv) #store as list obj.
## Check the following out:
# head(store_res) # list object with 2 slots
# str(store_res)
# store_res[[1]]

## Extract the values at column '10th percentile' at position nrep +1
## (= nrow command) to access avg value of all repetitions at the lowest row

zename<-'X10.percentile.training.presence.logistic.threshold' # check if the column name is the same !!!!!
# Achtung! 
# In MaxEnt 3.4.1 heisst Spalte 'X10.percentile.training.presence.Logistic.threshold'
zename

zecol <-which(colnames(store_res[[1]]) == zename) 
zecol

## a little head twister:
## access each element of your list, which is a whole dataset with store_res[[1]]
## and inside this list object = data.frame, access the element of
## the last row (nrow), and the column zecol] containing the threshold value
t_noBias   <- store_res[[1]][nrow(store_res[[1]]),zecol]
t_withBias <- store_res[[2]][nrow(store_res[[2]]),zecol]

t_noBias #the thresholds

t_withBias

all_bias <- rep(c(t_noBias,t_withBias), each=2)
all_bias

binary_thresh <- me_stack >= all_bias
plot(binary_thresh)

test05 <- me_stack >= 0.5 # check result with fixed threshold of default cut-off
plot(test05)

#testSeq <- me_stack >= seq(0.1, 0.7, length = 4)
#plot(testSeq)

# Extract relative probability values inside protected areas PAs
# check if overlay correct - only works when you have NOT aggregated the maps
compareRaster(PA_ras,me_stack[[1]]) 

ex <- which(values(PA_ras)==1)
# ex # gives Pointer to elements/ Index
ex_stack <- extract(x=me_stack, y=ex)
head(ex_stack)

boxplot(ex_stack,na.rm=T)

ex_binary <- extract(x=binary_thresh, y=ex)
head(ex_binary)

## In the following, we are summing the cells with 'TRUE'
ex_binary[1:2,]; colSums(ex_binary,na.rm=T) # 2 cmnds in line with ; !